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ABSTRACT 



Context. Three-dimensional (3D) radiative hydrodynamic model atmospheres of metal-poor late-type stars are characterized by cooler 
upper photospheric layers than their one-dimensional counterparts. This property of 3D model atmospheres can dramatically affect 
the determination of elemental abundances from temperature-sensitive spectral features, with profound consequences on galactic 
' chemical evolution studies. 

, Aims. We investigate whether the cool surface temperatures predicted by 3D model atmospheres of metal-poor stars can be ascribed 

to approximations in the treatment of scattering during the modelling phase. 

Methods. We use the Bifrost code to construct 3D model atmospheres of metal-poor stars and test three different ways to handle 
scattering in the radiative transfer equation. As a first approach, we solve iteratively the radiative transfer equation for the general 
^ , case of a source function with a coherent scattering term, treating scattering in a correct and consistent way. As a second approach, 

5_( ■ we solve the radiative transfer equation in local thermodynamic equilibrium approximation, neglecting altogether the contribution of 

' continuum scattering to extinction in the optically thin layers; this has been the default mode in our previous 3D modelling as well 

' as in present Stagger-Code models. As our third and final approach, we treat continuum scattering as pure absorption everywhere, 

which is the standard case in the 3D modelling by the CO^'BOLD collaboration. 

Results. For all simulations, we find that the second approach produces temperature structures with cool upper photospheric layers 
very similar to the case in which scattering is treated correctly. In contrast, treating scattering as pure absorption leads instead to 
^ ' significantly hotter and shallower temperature stratifications. The main differences in temperature structure between our published 

lO , models computed with the Stagger- and Bifrost codes and those generated with the CO^BOLD code can be traced to the different 

■ treatments of scattering. 

OA ' Conclusions. Neglecting the contribution of continuum scattering to extinction in optically thin layers provides a good approximation 

CO ' to the full, iterative solution of the radiative transfer equation in metal-poor stellar surface convection simulations, and at a much lower 

computational cost. Our results also demonstrate that the cool temperature stratifications predicted for metal-poor late-type stars by 
previous models by our collaboration are not an artifact of the approximated treatment of scattering. 

Key words. Hydrodynamics - Convection - Radiative transfer - Scattering - Stars: late-type - Stars: atmospheres - Methods: 
numerical 

1. Introduction of approx imate descriptions suc h as the mixing-length theory 

H . . (MLT) bv lB5hm -Vitensd d 19581) or the full spectru m of turbu- 

. , Model stellar atmospheres are indispensable tools for the quan- ign^e (FST) method bv lCanuto & Mazzitellil JMl, which are 
titative mterpretation of observed stellar spectra and colours, ^jj hampered by a number of free parameters. 
Synthetic stellar fluxes from model stellar atmospheres can for 

instance be compared with observations to infer the physical Recent years, on the other hand, have seen a rapid devel- 
properties and chemical compositions of distant stars. opment of numerical three-dimensional (3D) radiative hydrody- 

The vast majority of theoretical model atmospheres of solai-- namic simulations of stellar surface convection. In these simula- 
and late-type stai's in use today ai-e still computed under the tions, the gas flows in the highly stratified outer layers of stars are 
assumption of one-dimensional (ID) geometry, flux constancy, modelled by solving the hydrodynamics equations of mass, mo- 
and hydrostatic equilibrium. Yet, convective energy transport - mentum, and energy conservation and accounting for the energy 
which in late-type stars significandy affects the photospheric exchanges between matter and radiation via radiative transfer In 
temperature stratification - can be modelled in ID only by means this framework, convective motions arise and self-organize nat- 

urally without the need to introduce adjustable parameters. One 

Send offprint requests to: R. Collet, of the main goals of 3D modelling of stellar surface convection 

e-mail: remo@mpa-garching.mpg.de is therefore to provide a more realistic description of the phys- 
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ical structure o f late-type stellar atmospheres. Since the earl y 
studies by, e.g., 'Nordlund' ('1982|), iNordlu nd & DravinsI (Il990l) . 
and Steffen et al. (1989), convection simulations have become 
increasingly more sophisticated not only from the point of view 
of the implemented numerical algorithms, but also in terms of 
input physics, allowing for quantitative comparisons with obser- 
vations. Surface convection simulations have been successfully 
applied to studying the properties of granulation at the surface 
of the Sun (e.g. jS tein & Nordlund 1998; Carlsson et al. 2004; 
Danilovic et al.'2008; Pereira et al. 2009b; Nordlund et al..2009i: 
Wed emever-Bohm & Rouppe van der Voort 2009 ) and other 
late- type stars (e.g. Allen de Prieto et alj 120021 ; IRamfrez et al.l 
5009 ; Chiavassa et al. 20ldr Theoretical predictions of pho- 
tospheric temperature stratifications from 3D model stellar at- 
mospheres have been tested in a number of cases against 
observed centre-to-limb variations, showi ng in general excel- 
lent agreem ent with the measurements ( Aufdenberg et al."2005'; 



Bigot et al 



Chiavassa et al 



2006,; |Kpesterke et al., 
20 iK 



■2008; ,Pereii-a et al.. ,2009ai: 



Though not yet commonplace, 3D model atmospheres 
have also started to be employed in spectroscopic analyses 
for the determination of stellar elemental abundances (e.g . 



Asplund et al.l 1999', l2003l, 12009"; 'Asplund & Garcia Per e J20oT 
Collet et al. 2006, 2007; Caffau et al. 2010; Behara et al.lboiC: 



Gonzalez Hernandez et al.l l2010). The temperature and density 



inhomogeneities and the velocity fields present in 3D hydrody- 
namic simulations but absent in ID hydrostatic models are in 
general responsible for differences between 3D and ID analyses 
in terms of predicted shapes and strengths of synthetic spectral 
lines and in terms of abundances inferred from observed stel- 
lar spectra. For very metal-poor stars, in addition, another aspect 
contributes to differences between 3D and ID spetroscopic anal- 
yses. Three-dimensional surface convection simulations predict 
in fact much cooler temperature stratifications with steeper gra- 
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dients than ID hydrostatic models in the upper layers of ver 
metal-poor late-ty pe stellar atmospheres jAsplund et al.l 
ICollet et al.ll2007h . In the upper layers of 3D hydrodynamic sim- 
ulations, the temperature stratification is determined by the com- 
petition between radiative heating and adiabatic cooling owing 
to the expansion of gas overshooting into the atmosphere from 
the top of the convection zone. In solar-metallicity 3D simula- 
tions, spectral lines contribute enough heating in the upper at- 
mosphere to maintain the temperature stratification close to a 
purely radiative equilibrium one, as in ID models. At very low 
metallicities, on the contrary, the scarcity and overall weakness 
of line opacity sources and the cooling effect owing to the ex- 
pansion of gas above granules act together to keep the average 
temperature low in the upper photosphere. In corresponding sta- 
tionary ID hydrostatic models, the cooling component associ- 
ated with gas flows is missing and balance is reached at higher 
temperatures. From the point of view of spectroscopic analy- 
ses, these structural differences between 3D and ID models can 
lead to large differences (up to about 1 dex) in terms of derived 
abun dances wheneve r temperature-sensitive spectral features are 
used (lAst?lundll200"5h . While there is general qualitative agree- 
ment on the above results, the predicted average equilibrium 
temperature stratifications of 3D model metal-poor atmospheres 
from different groups often vary significantly depending on the 
adopted codes and input physics. These differences are partly re- 
sponsible for discrepancies in terms of computed 3D- IDcorrec- 
tions to elemental abundances from different studies (Bo nifaciol 
[2010; Dobrovolskas et al. 2010; Ivanauskas et al. 2010). 

We propose that different treatments of scattering in the solu- 
tion of the radiative transfer equation are responsible for the dif- 



ferences among various convection codes in terms of predicted 
photospheric temperature stratifications. We explore for the first 
time the effect of coherent isotropic continuum scattering on the 
energy balance in the upper layers of 3D stellar surface convec- 
tion simulations of metal-poor late-type stars. We compute simu- 
lations where we iteratively solve the radiative transfer problem 
for a source function that includes a coherent scattering term, 
and compare them with other simulations where scattering is 
treated by means of two commonly adopted approximations, by 
either neglecting it altogether in the optically thin layers or by 
including it as true absorption everywhere. 



2. Methods 

2.1. Surface convection simulations 

We carried out surface convection simulations of two metal-poor 
red giants and one metal-poor turnoff star using t he 3D, time- 
dependent, radiation-hydro dynamic Bifrost code dHavek et al.l 
2010; Carlsson et al. 2010; Gudiksen et al., in prep.). The mass, 
momentum, and energy conservation equations are discretized 
using a high-order finite-difference scheme (6* order deriva- 
tives, 5* order interpolations) and solved on a staggered Eulerian 
mesh for a representative, 3D, rectangular, volume located 
across the optical surface. We employ a numerical resolution of 
240^ for the mesh, with five ghost zone layers at the top and 
the bottom of the simulation box. As for the boundary condi- 
tions, we adopt periodic boundary conditions horizontally and 
open boundaries vertically. At the bottom of the simulation do- 
main, the inflowing material is set to have constant entropy and 
pressure, while the outflowing material is free to carry entropy 
fluctuations out. In the last physical layer at the top of the simula- 
tion box, the internal energy per unit mass is forced not to deviate 
by more than 10% from its average value, for stability reasons. 
Also, for inflowing gas, the internal energy per unit mass is then 
set to its minimum value from the uppermost five physical layers 
at the top. 

The stellar parameters of the simulations are summarized in 
Table [1] For each si mulation, we assume a solar composition 
(Asplun d et al.ll2005h with the abundances of all metals scaled 
proportionally to the iron abundance. We chose the physical 
dimensions of the simulations to be sufficiently large to allow 
about ten major granules to develop at any one time in the do- 
main and to span about twelve pressure scale heights vertically. 
In terms of Rosseland optical depth, the simulations cover the 
range -4SlogTRoss^6. As initial conditions, we used snapshots 
from previous simulation sequences by Collet et al. (2007) and 
lAsplund & Garcia Perezi (120011) and converted them to the de- 
sired numerical resolution. We then let the simulations achieve 
complete thermal relaxation by running them over a period of 
one or two convective overturn timesQCiven that the fundamen- 
tal stellar parameters of the simulations were not varied and that 
the changes owing to updates of the input micro-physics are rela- 
tively small, relaxation is rapidly achieved during this time span. 



2.2. Radiative transfer 

An accurate representation of radiative energy exchange is es- 
sential for a realistic description of the temperature stratification 



' We define here the simulation's convective overturn time as the 
characteristic time it takes for a fluid parcel to travel through the en- 
tire depth of the simulation domain, from the bottom to the surface and 
then back. 
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in stellar surface layers. We describe the interaction between the 
plasma and radiation via a radiative heating term in the energy 
conservation equation: 

Grad^ r r XA(lA-SA)dAdn, (1) 

JnJo 

where xa is the total monochromatic extinction coefficient (per 
unit path length), sum of the absorption k,i and scattering cta 
extinction coefficients, Sa and are the monochromatic source 
function and intensity, respectively, and the integral is performed 
over the whole wavelength range and over the whole solid angle. 
Under the assumption of an isotropic source function and extinc- 
tion coefficients and of a static medium, integral ^ is equivalent 
to 

erad = 47r I XA(JA-SA)dA, (2) 

Jo 

where Ja - j;; J^h dO. is the monochromatic mean inten sity. 

We adopt a realistic equation of state (iMihalas et al.lll988l) . 
which accounts for the effects of excitation, ionization, and dis- 
sociation of the 15 most abundant elements and of H2 and 
H^, continuous opacities from Gus tafs son et al.l (Il975h and 
Trampedach (2011, in prep.), and sampled line opac ities from 
B. Plez (priv. comm.; see also iGustafsson et al.ll2008|). The full 
list of references for the opacity sources is given in lHavek et aP 
(Hop). In our simulations, we compute the radiative heating 
rates by solving the radiative transfer equation at each time step 
and at each grid-point in three-dimensions on short character- 
istics along 24 rays. The integral over solid angle in Eq. [T] is 
then appr oximated by a weighted Carlson's A4 quadrature sum 
(ICarlsonr i963) of the Ii — Sa terms at these specific directions. 
The implemented radiative transfer method dHavek et al..,2010.) 
allows for a solution of the radiative transfer equation with co- 
herent scattering by means of an iterative Gauss-Seidel acceler- 
ation scheme (Truiillo Bueno & Fabiani Bendicho 1995). In the 
latter case, the monochromatic source function is expressed in 
the form Sa - £aBa + (l-eO-^i^ where £a=ka/xa is the photon de- 
struction probability, Ba is the Planck function, and eaBa is the 
thermal emissivity. Below, we will refer to the solution of the 
radiative transfer equation for a source function with a coherent 
scattering term as the standard case. With a source function of 
the above form, it is easy to show that the radiative heating rate 
(|2|l can be expressed as 

emd=47r KA(JA-BA)dA. (3) 

Jo 

Equation (O emphasizes that for radiative transfer with coher- 
ent scattering only absorption processes effectively contribute to 
radiative heat exchange. The consequences of this will be dis- 
cussed in more detail below, in Sect. |4] 

Continuum scattering plays a negligible role in shaping 
the photospheric stratification of solar-type stars (Havek et al.l 
I2OIO1) . In the much thinner atmospheres of metal-poor giants, 
however, where the electron number density is also lower, the 
absorption opacities as well as the rates of collisional processes 
that can thermalize the radiation field are strongly reduced and 
scattering becomes an important opacity source. As an exam- 
ple, Fig.[T]illustrates the dependence of the destruction probabil- 
ity €a on wavelength and depth in the mean stratification taken 
from a m etal-poor red giant surface convection simulation by 
ICollet et a l. (2007). Deep below the optical surface (located at 
» Mm in the figure), absorption processes dominate the inter- 
action between gas and photons, and the radiation field is fully 



Table 1. Parameters of the 3D hydrodynamic simulations. 



[K] 


[cgs] 


[Fe/H] 


x,);,z-dimensions 
[Mm] 


Time^ 
[hrs]~ 


[hrs] 


5140 


2.2 


-3.0 


1150x1150x430 


28 


20 


5080 


2.2 


-2.0 


1150x1150x430 


28 


20 


6510 


4.04 


-3.0 


21.4x21.4x8.35 


1.3 


0.3 



3.5 
1.5 
0.15 



" Temporal average of the emergent effective temperatures. 
* Time span of the simulations. 
'' Average granule lifetime. 

'' Radiative heating timescale when switching to scattering as true 
absorption in the upper photosphere. 

Table 2. Bin boundaries: a wavelength A is assigned to a certain 
opacity bin depending on the Rosseland optical depth at which 
the monochromatic optical depth ta is equal to one. The table 
shows the Rosseland optical depth ranges spanned by individual 
bins in the present study. 



Bin 


logTRoss(T,, 


= 1) range 


1 


-Hoo 


-0.5 


2 


-0.5 . . . 


-1.5 


3 


-1.5 ... 


-2.5 


4 


-2.5 ... 


— CO 



thermalized (ei=l). Above it, continuum scattering - in partic- 
ular, Rayleigh scattering by H i in the UV and blue part of the 
visible spectrum and electron scattering at longer wavelengths 
- becomes increasingly important as an opacity source, and the 
destruction probability £a is closer to zero for a large part of the 
spectrum. 

2.3. Opacity binning 

Solving the radiative transfer equation for a large number of 
wavelengths in 3D convection simulations is a computation- 
ally demanding task. We therefore use here the opacity binning 
or multigroup method (Nordlund 1982) to approximate the de- 
pendence of opacities on wavelength. The problem of resolv- 
ing the more than 100000 wavelengths of the original opacity 
sampling data is reduced to computing the radiative transfer for 
only a small number of mean opacities. Monochromatic contin- 
uous plus line opacities are sorted into four representative bins 
(groups) according to their strength. Bin membership is deter- 
mined by the height where monochromatic optical depth unity 
is reached on the simulations' mean temperature-density strat- 
ifications. The latter are computed beforehand from thermally 
relaxed simulation sequences^ As the reference height scale we 
chose here the Rosseland optical depth computed for each sim- 
ulation's mean stratification. The bin boundaries for the three 
simulations considered here are shown in Table |2] 

We followed the formalism of Skartlien (2000) to compute 
the mean opacities and integrated thermal emissivities for each 
bin. In computing the mean opacities, we distinguished between 
diffusion and streaming regimes, depending on whether the av- 
erages are computed in the optically thick or in the optically thin 
layers of the mean temperature-density stratifications, respec- 

' More specifically, we compute here the mean stratifications by av- 
eraging the simulation snapshots over time and over surfaces of constant 
column mass density. 
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Wavelength [A] 

Fig. 1. Photon destruction probability as a function of wavelength and geometrical depth in the mean stratification of a metal-poor 
([Fe/H]=-3) red giant surface convection simulation. Darker shades indicate lower photon destruction probability. The destruc- 
tion probabilities are here computed assuming that only continuous processes contribute to scattering extinction. Spectral lines 
(appearing as white vertical stripes in the upper part of the figure) are treated in pure absorption. 



lively. In the diffusion limit, we define the extinction coefficient 
in bin / as the Rosseland group mean 



(4) 



where AAi denotes the set of wavelength intervals assigned to bin 
i. In the streaming limit on the other hand, we approximate the 
extinction coefficient in bin / with an intensity-weighted mean 



JaAi I JaAi 



(5) 



where is the monochromatic mean intensity computed for 
the mean (plane-parallel) atmospheric temperature-density strat- 
ification and averaged over the whole solid angle. In general, 
we define the mean extinction coefficient as a height-dependent 
linear combination of the values in the diffusion and streaming 
limits by using an exponential bridging 

X^^e-f^X^ + (l-e-'^")xf (6) 

where Tf is the Rosseland group mean optical depth computed 
for the mean stratification and / is a constant factor that we sim- 
ply set equal to 30 in the present study. The exponential bridg- 
ing allows for a smooth transition from streaming to diffusion 
regime, while the choice of / = 30 ensures that the stream- 
ing regime component is adequately suppressed below the opti- 
cal surfaceQ Bin opacities for arbitrary temperature and density 
pairs are then simply determined by constant extrapolation of /pp 
on surfaces of constant Rosseland optical depth of the opacities 
calculated for the mean stratification. 

In order to investigate the effective role of scattering on the 
radiative heating balance in 3D surface convection simulations, 
we considered three methods to handle it in the radiative transfer 
calculations. 

As a default, we solve the radiative transfer equation consis- 
tently for the general case of a source function with a coherent 
isotropic continuum scattering term. We will refer to this as the 
standard case. More precisely, for each bin we assume a source 
function of the form 



S i ^ (£B)i + (I - edJi, 



(7) 



^ We also performed test calculations with opacity tables computed 
for a value of f=2 of the bridging parameter and found that the resulting 
stratifications from the simulations are essentially unaltered. 



where 7, is the mean intensity for bin /, 



(eB),- = f 



kaBa dA 



is the integrated thermal emissivity, 
Kf ^ f KAJ^dAl f f^dA 

JaA, I JAAi 



and 



erf = r CTAJ'/dAl f J^UA 

JAA, ' JAA 



'A A 
)AAi I JAA; 



(8) 



(9) 



(10) 



are the streaming-regime mean absorption and mean scattering 
extinction coefficients, respectively, and 



(11) 



is the mean photon-destruction probability. The mean intensity 
Ji is determined by solving the radiative transfer equation for in- 
te nsity at each time-st ep with the iterative method implemented 
bv lHaveketaD(l20Toh and averaging over solid angle: 



(12) 



where 7,^ is the intensity in bin / and along direction j, and the 
<jj/s are the weights of the Carlson's quadrature. Finally, with 
the opacity binning approximation the radiative heating rate Q 
becomes 



firad =^^A'^w>(/,7 



•5,■) 



:47r2 XiiJi-Si). 



(13) 



As a second approach, we approximated the monochromatic 
source function with the Planck function {Sa - B x) through- 
out the simulation domain and completely excluded the con- 
tribution of continuum scattering from the calculation of the 
mean extinction coefficients in the streaming regime in each 
bin. This is essentially the same approxima t ion of scattering 
as the one adopted by iNordlund & DravinsI (1 19901) (see also 
Skaxtlien 2000, Sect. 5.2.1) and subsequently in stellar sur- 
face convection sir nulations carried out by our group (e.g., 
lAsplund et aT]ll999HAsplund & Garcia Peredl200lLrCollet et all 
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l2007h with the lStein & Nordlundl (IT998h code, or bvlCoUet et alJ 
(l2009h with the Stagger-Code (iNordlund & Gals gaardlll995r 
Compared with the standard case, this method is faster and com- 
putationally less demanding. 

With this approximation, the integrated source function in 
bin i is given by 



Si = Bi 



j 

Jaa, 



BidA 



(14) 



and the mean extinction coefficient in the streaming regime by 

;rf*=^/= r K.jYdAl f J^dA. (15) 
Jaa, I Jaa, 

In the opacity binning approximation, the radiative heating rate 
(H]) becomes 



(16) 



where 7, is again the mean intensity for bin /. With our assump- 
tion ( fTsT i for the extinction coefficient in the streaming regime, 
the radiative heating rate in the optically thin layers of the con- 
vection simulations approaches 



.47r^ KiiJi-Bi). 



(17) 



Finally, as a third approach, we still assumed a Planckian 
source function, but this time included the contribution of scat- 
tering to the total extinction in both the diffusion and the stream- 
ing regimes, effectively treating scattering as true absorption ev- 
erywhere in the simulation domain. This corresponds to the ap- 
proximation adopted in stellar surface convection simulations by 
the CO^BOLD collaboration (H. G. Ludwig, priv. comm.). 

For each set of stellar parameters indicated in Table [1] start- 
ing from the same initial snapshot, we computed three simula- 
tion sequences - one for each of the different treatments of the 
radiative transfer considered above - and followed the evolution 
of the temperature stratification with time in all three cases. 

3. Results 

Figure|2]compares the resulting temperature structures from the 
[Fe/H]=-3 red giant surface convection simulations computed 
for the three cases: (i) self-consistent solution of the radiative 
transfer solution for a source function with a coherent scatter- 
ing term (standard case); (ii) Planckian source function and no 
scattering contribution to extinction in optically thin layers; (iii) 
scattering treated as true absorption everywhere. The tempera- 
ture stratification resulting from the assumption of a Planckian 
source function and of no contribution of scattering to extinc- 
tion in optically thin layers is very similar to the one predicted 
in the standard case (Fig. |2l left panel). In particular, the over- 
all temperature distrib ution with op tical depth is consistent with 
the previous results bv lCoUet et al.l (12007.) . As shown in the left 
panel of Fig. [3l the snapshot-by-snapshot temperature diff'er- 
ence between the mean stratifications in the standard case and in 
the no-scattering-in-streaming-regime case is typically less than 
50 K at all depths. Moreover, not only the mean stratifications 
resulting from these two treatments of scattering in the solution 
of the radiative transfer equation are virtually identical, but the 
two simulations also maintain a good correlation in terms of 3D 
temperature and density structures as well as velocity fields af- 
ter several (^ 12) hours of simulated stellar time (the average 
granule lifetime is about 20 hours). 



Conversely, including scattering as true absorption while as- 
suming a Planckian source function causes a rapid heating of the 
upper photospheric layers immediately after starting the simula- 
tion (Fig. [3] right panel). With respect to the standard simula- 
tion, the temperature at logTRo,ss=-4 rises by 500 K in about 
six hours of simulated stellar time. As the simulation evolves, 
the absolute temperature difference with respect to the standard 
case settles to about 600 K at that depth (Fig. |2] right panel). 
The temperature rise with respect to the standard case is approx- 
imately proportional to 1 - exp(-f/frad), where f^d is the char- 
acteristic radiative heating timescale of the upper photospheric 
layers when scattering is switched to true absorption. In partic- 
ular, for the [Fe/H]=-3 red giant simulation, frad~3 hours; the 
characteristic radiative heating timescales for the other simula- 
tions are listed in Table [1] Mixing caused by bulk gas flows has 
an effect on the temperature adjustments in the atmospheres, but 
on longer timescales. At greater depths the differences are less 
pronounced and amount to about 300 K at logTRoss=-3 when 
the temperature balance is achieved and to less than 50 K be- 
low logTRoss=-2. This dependence of the heating on depth also 
results in a shallowing of the vertical temperature gradient in 
the upper photosphere. Near the optical surface, the temperature 
differences are negligible and the effective temperature of the 
simulation is practically unaltered. 

Figure |4l left panel, shows the temperature distribution at 
three different optical depths for the relaxed [Fe/H]= -3 red 
giant simulations computed with the three different treatments 
of scattering. Again, the distribution returned by the simulation 
with proper treatment of coherent scattering agrees very well at 
all optical depths with the one predicted by the standard case, but 
differs significantly from the one predicted in the scattering-as- 
absorption case. In particular, at an optical depth of log rRo.ss=-3, 
where the differences between the scattering-as-absorption and 
coherent scattering cases are more evident, the latter predicts 
a skewed temperature distribution peaking at 3400 K, with a 
longer tail at higher temperatures and with the hint of a smaller 
secondary peak at 4000 K; the simulation that assumes scattering 
in absorption instead yields a distribution skewed towards higher 
temperatures and peaking at 4200 K and with a less pronounced 
peak at around 3400 K. The overall shape of the distribution is 
determined by the competition between dynamic and radiative 
heating/cooling processes. Radiative heating/cooling reduces the 
skewness of the distribution by attenuating temperature fluctua- 
tions. The secondary peak as well as the low-temperature tail 
of the distribution in these optically thin layers are mainly con- 
tributed by cool inflowing gas from the top of the simulation 
domain. At greater depths, the differences in temperature distri- 
bution between the three simulations become less and less sig- 
nificant. 

Figure|5]shows the resulting distribution of temperature with 
optical depth for the metal-poor ([Fe/H]=-3) turnoff-star sim- 
ulations computed in the standard (left panel) and scattering- 
as-absorption (right panel) cases. As for the metal-poor red gi- 
ant simulations at [Fe/H]=-3, the stratification predicted by the 
standard case agrees well with the one predicted by the simu- 
lation that assumes a Planckian source function and no scatter- 
ing contribution to extinction in optically thin layers, while the 
case where scattering is treated as pure absorption everywhere 
produces a hotter and shallower temperature stratification in the 
upper photosphere. The corresponding histograms in the central 
panel of Fig. |4] illustrate in more detail the temperature distri- 
butions at various depths in the three cases. Again, the distri- 
butions in the standard and in the no-scattering-in-streaming- 
regime cases are very similar at all depths, while the tempera- 
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Fig. 2. Grey shaded areas: distribution of temperature as a function of Rosseland optical depth in two representative snapshots of 
two surface convection simulations of a metal-poor red giant at [Fe/H]=-3. The snapshots were taken 28 hours after the start of the 
simulations. Darker shades indicate temperature values with higher frequency of occurrence. Left panel: temperature distribution 
from the simulation with self-consistent solution of the radiative transfer equation for a source function with a coherent scattering 
term (standard case). Continuous blue line: mean temperature stratification from the corresponding simulation assuming a Planckian 
source function and neglecting the contribution of scattering in optically thin layers, averaged over time and surfaces of constant 
Rosseland optical depth. Right panel: temperature distribution for the radiative transfer solution assuming a Planckian source func- 
tion and treating scattering as true absorption everywhere. Continuous red line: corresponding mean temperature stratification. 
Dashed lines (both panels): mean temperature stratification for the standard-case simulation. 
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Fig. 3. Temporal evolution of the mean temperature stratification from the two metal-poor red giant ([Fe/H]=-3) surface convection 
simulations considered in Fig. |2] The curves show the differences of the mean temperature with respect to the standard case (self- 
consistent solution of the radiative transfer equation for a source function with a coherent scattering term) as a function of time. The 
instantaneous mean temperature stratifications are computed on surfaces of constant Rosseland optical depth. The total temporal 
coverage is 28 hours, with snapshots taken every 10 minutes. Darker curves indicate later snapshots in the two sequences. Left panel: 
temporal evolution of the temperature differences for the simulation computed assuming a Planckian source function and neglecting 
the contribution of scattering to total extinction in optically thin layers. Right panel: corresponding evolution for the simulation 
computed under the assumption of Planckian source function with scattering treated as true absorption everywhere. 



ture distribution for the simulation in which scattering is treated 
as pure absorption clearly departs from the other two in the up- 
permost photospheric layers. 

The temporal evolution of the mean temperature stratifica- 
tions is shown in Fig.|6] Because of the stronger surface gravity 
and smaller spatial dimensions, the adjustment timescales are 
shorter for the turnoff-star simulations than for the red giant se- 
ries (f,ad~10 minutes for the turnoff sequence). The outer layers 
of the simulation in which scattering is treated as pure absorption 
heat up rapidly, and reach a stable configuration within approxi- 



mately a few tens of minutes of stellar time, with the temperature 
being on average 600 K and 200 K higher than in the simula- 
tion with coherent scattering at log tross=-4 and at log tross=-3, 
respectively. Because the dynamical adjustment timescales are 
short, the temperature and density perturbations arising from the 
different treatments of scattering in the radiative transfer grow 
rapidly and the correlation between simulations from the two 
different cases degrades relatively fast. This is partly the rea- 
son why the snapshot-by-snapshot temperature differences be- 
tween the mean stratifications in the standard and scattering- 
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Fig. 4. Histograms of temperature in the upper photospheric layers of the surface convection simulations at three different Rosseland 
optical depths and for the three different treatments of scattering in the radiative transfer solution. Left panel: results for the red giant 
at [Fe/H]=-3; central panel: turnoff star at [Fe/H]=-3; right panel: red giant at [Fe/H]=-2. The histograms are constructed from 
snapshots taken at the end of the simulation sequences, for the coherent scattering case {grey bars), Planckian source function and no 
scattering extinction in optically thin layers {blue curves), and Planckian source function with scattering treated as true absorption 
{red curves). 



in-absorption cases above the optical surface are slightly larger 
than for the [Fe/H]=-3 red giant simulations, although the dif- 
ferences typically remain below 100 K. 

We also carried out a series of simulations of a metal-poor 
red giant at slightly higher metallicity ([Fe/H]=-2) for the three 
different treatments of scattering. The results (Fig. [T] and [8l 
Fig. m right panel) are qualitatively similar to the [Fe/H]=-3 
simulations. The scattering-in-absorption case also produces a 
hot stratification in the uppermost part of the photosphere, but 
the deviations from the standard case are smaller (a!250 K at 
log tro,ss=-4). Excluding scattering from the total extinction in 
the optically thin layers again leads to a stratification that agrees 
well with the coherent scattering case at all depths. However, 
in the upper part of the photosphere, the stratification from the 
simulation assuming a Planckian source function and no contri- 
bution of scattering to extinction in optically thin layers is on 
average slightly cooler (by 100 K at logTRos.s=-4) than the one 
predicted in the standard case. 

4. Discussion 

Our results indicate that assuming a Planckian source function 
and excluding scattering in optically thin layers is a good ap- 
proximation for estimating radiative heating rates in metal-poor 
late-type stellar surface convection simulations. At first sight, 
this is an intuitive result to some degree, since coherent scatter- 
ing does not contribute directly to the energy balance in the up- 
per photosphere: its primary effect is, in fact, to redirect photons 



to other directions without changing their energy. On the other 
hand, when scattering processes coexist with absorption pro- 
cesses in stellar atmospheres (that is, when scattering and pure 
absorption opacities are comparable), they affect the form of the 
source function to an appreciable degree and, consequently, play 
a crucial role in determining the actual magnitude of mean in- 
tensities and of radiative heating rates. 

The non-local nature of scattering and the presence of tem- 
perature and density inhomogeneities in late-type stellar photo- 
spheres renders the problem of determining the effective feed- 
back of scattering on the energy balance very complicated. It is 
therefore necessary to carry out numerical simulations to verify 
our hypothesis that excluding the contribution of scattering to 
extinction in the optically thin layers may work as a good ap- 
proximation to the correct solution. With hindsight, we can use 
our results to justify why our approximation works. 

Below the optical surface, the assumption of a Planckian 
source function is relatively accurate because of the dominance 
of absorption that contributes to the thermalization of the ra- 
diation field. Proceeding upwards, the opacity changes from 
absorption-dominated to scattering-dominated. In metal-poor 
late-type stars (and particularly in giants), this transition is fairly 
rapid, as we have seen (Fig. [T]i. In the optically thin regions of 
the photosphere, where scattering processes dominate, the mean 
intensity is to first order set by the physical conditions in the lay- 
ers located underneath, near the optical surface. By neglecting 
the contribution of scattering to the opacity in surface layers and 
assuming that the sole role of scattering there is to redirect pho- 
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Fig. 5. Left panel: temperature distribution with Rosseland optical depth from a representative snapshot of the [Fe/H] = -3 turnoff-star 
surface convection simulation computed solving the radiative transfer equation self-consistently for a source function with a coherent 
scattering term (standard case). Darker shades indicate higher probability. Blue curve: mean stratification from the corresponding 
simulation computed assuming a Planckian source function and no contribution of scattering to total extinction in optically thin 
layers. Right panel: corresponding temperature distribution with Rosseland optical depth from the simulation assuming a Planckian 
source function and scattering as true absorption everywhere (mean temperature stratification shown in red). Dashed lines (both 
panels): mean temperature stratification from the standard-case simulation. All snapshots were taken 80 minutes after the start of 
the simulations. 
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Fig. 6. Temporal evolution of the mean temperature stratification from the two metal-poor ([Fe/H]=-3) turnoff-star surface con- 
vection simulations considered in Fig. |5] The curves show the differences of the mean temperature with respect to the standard 
case (self-consistent solution of the radiative transfer equation for a source function with a coherent scattering term) as a function 
of time. The total temporal coverage is 80 minutes, with snapshots taken every minute. Darker curves indicate later snapshots in 
the two sequences. Left panel: temporal evolution of the temperature differences for the simulation assuming a Planckian source 
function and no contribution of scattering to total extinction in optically thin layers. Right panel: corresponding evolution under the 
assumption of Planckian source function with scattering treated as true absorption. 



tons, we can provide a rough estimate of the mean intensity in the 
upper photospheric layers and simulate these non-local effects in 
the radiative transfer even though we have adopted a Planckian 
source function everywhere for the calculations. Naturally, the 
mean intensities predicted by the coherent scattering case and 
by the approximated treatment with no scattering contribution 
to extinction in the streaming regime disagree to some extent. 
Small differences also arise in practice from the different weight- 
ing of opacities in the calculation of the mean extinction coeffi- 
cients xi ™d xi* (Eq- 12] and [TSl l and can contribute to small 
deviations in terms of predicted radiative heating rates. 



Let us now consider what happens when we start from a 
snapshot of the simulation assuming a Planckian source function 
and no scattering contribution to extinction in optically thin lay- 
ers and switch to the scattering-as-absorption case. Changes to 
the total radiative heating rate in the uppermost layers are most 
significant in the first bin that groups the wavelengths for which 
the formation region is located near the optical surface, while 
they are largely negligible for the remaining three bins. This is il- 
lustrated in Fig.|9l which shows the comparison between the hor- 
izontally averaged instantaneous radiative heating rates in indi- 
vidual bins computed for a representative snapshot of the metal- 
poor red giant simulation at [Fe/H]=-3 that assumes no scat- 
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Fig. 7. Temperature distribution with Rosseland optical depth at the surface of a metal-poor red giant at [Fe/H]=-2. Left panel: pre- 
dicted distribution from a representative snapshot of the simulation computed solving the radiative transfer equation self-consistently 
for a source function with a coherent scattering term (standard case). Blue curve: mean stratification from the corresponding sim- 
ulation computed assuming a Planckian source function and excluding scattering from extinction in optically thin layers. Right 
panel: corresponding distribution from a simulation assuming a Planckian source function with scattering as true absorption (mean 
stratification shown in red). Dashed lines (both panels): mean stratification from the standard-case simulation. All snapshots were 
taken 28 hours after the start of the simulations. See Fig.|2]for more details. 



800 



<u 600 
u 



^ 400 

Q 
0) 

3 200 
o 

(U 
Q. 

E 
<u 



T^„=5080K logg = 2.2 [Fe/H] = -2 
S=B, x''=f'' - S=eBH-(1-e)J 



-4 



-2 
log -^Ross 



800 



<u 600 
o 



^ 400 

Q 
(U 

3 200 
o 

0) 
CL 



T^„=5080K logg = 2.2 [Fe/H] = -2 
S=Q, x^=K^+a'' - S=eBH-(1-e)J 




-4 



-2 
log -^Ross 



Fig. 8. Temporal evolution of the mean temperature stratification in the two metal-poor red giant simulations at [Fe/H]=-2 from 
Fig. Q The curves show the differences of the mean temperature with respect to the standard case (self-consistent solution of the 
radiative transfer equation for a source function with a coherent scattering term) as a function of time. The total temporal coverage is 
28 hours, with snapshots taken every hour. Left panel: temporal evolution for the simulation assuming a Planckian source function 
and no contribution of scattering to the total extinction in optically thin layers. Right panel: corresponding evolution under the 
assumption of Planckian source function with scattering treated as true absorption. See Fig.|3]for more details. 



tering contribution to extinction in the upper layers, before and 
after switching to the scattering-as-absorption case. Where the 
contribution of scattering to extinction in the streaming regime 
is neglected, the radiative heating rate in the first bin and in op- 
tically thin layers can be approximated by 

QUi«^7^k{{Ji-Bi). (18) 

The steepness of the temperature gradient ri ear the optical sur- 
face (lAsplund et al.lll999[ ICollet et al.ll2007h in the simulations 
where the contribution of scattering to extinction in optically thin 
layers is neglected ensures that the 7i - Bi split is generally pos- 
itive in the upper photosphere. As we will discuss below, this 
leads to radiative re-heating that maintains the temperature away 
from a pure adiabatic stratification. When scattering is switched 



to pure absorption, the opacity in the first bin significantly in- 
creases with regard to other bins, particularly in the upper pho- 
tospheric layers. The value of the mean intensity Ji is controlled 
by the conditions around the optical surface, where the changes 
in terms of opacity are small, and remains roughly the same, so 
that the Ji - B\ split is also unchanged to first approximation. 
The radiative heating rate in the first bin therefore becomes 

e;;d,i - ^^i>^x + <)iJi - > qUi > o, (i9) 

which explains why the temperature increases in the upper pho- 
tosphere when scattering is treated as pure absorption. 

Haveketal. (2010) have shown that continuous scattering 
plays a negligible role in shaping the temperature stratification in 
simulations of solar-type, solar-metallicity, stars. However, the 
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Fig. 9. Blue curves: horizontally averaged radiative heating rates for individual bins as a function of geometrical depth in a snapshot 
of the [Fe/H]=-3 red giant simulation computed assuming a Planckian source function and excluding the contribution of scattering 
to total extinction in optically thin layers. For clarity, only the relevant (upper) part of the simulation is shown; radiative heating 
rates at the cooling peak are unchanged when switching between scattering treatments. Red curves: horizontally averaged heating 
rates for the same snapshot immediately after switching to scattering as true absorption. 




Fig. 10. Grey shaded area: distribution of temperature versus gas pressure in two representative snapshots of metal-poor red giant 
surface convection simulations at [Fe/H]=-3 {left panel) and [Fe/H]=-2 {right panel). The simulations were computed with the 
radiative transfer solver for the coherent scattering case. Darker shades indicate higher frequency of occurrence. Dashed lines: mean 
stratifications, averaged over time and surfaces of constant Rosseland optical depth. Thin grey lines: contours of constant entropy. 



situation considered here is qualitatively different from theirs: in 
the upper photospheres of stars similar to the Sun, continuum 
scattering represents only a negligible fraction of the total opac- 
ity, so that excluding it or including it as true absorption in opti- 
cally thin layers does not lead to significant differences in terms 
of the predicted temperature and density stratifications. 

Throughout our analysis, we treated li ne opacities as pure 
absorption when computing the opacity bins. lHavek et alJ (1201 Ol) 
have also considered the case where line scattering is treated as 
such in the radiative transfer calculation s, estimating line photon 
destruction probabilities by means of the' van Regem orte^(ll962h 
approximation. According to their calculations, the overall im- 
pact on the temperature stratification is fairly small, with differ- 
ences of less than 40 K for log tross ^ -4 with respect to the 



case where line scattering is treated as pure absorption. It would 
be interesting to examine whether the effects of line scattering 
on heating rates are more significant for metal-poor stellar sur- 
face convection simulations, but this goes beyond the scope of 
the present discussion; we therefore reserve the study of these 
effects for future work. 

Another limitation to the accuracy of the calculations is the 
opacity binning approximation of non-grey radiative transfer. 
For the current analysis, we used a somewhat crude choice for 
the binning with only four opacity groups and no explicit divi- 
sion in the sense of wavelength. Allowing for more than four 
bins and grouping opacities according not only to strength but 
also to wavelength would increase the accuracy of the calcula- 
tions, but it is unlikely that it would affect the main conclusions; 
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the exact value of the difference between the temperature strat- 
ifications predicted with various treatments of scattering would 
depend on the details of the binning, but the deviations from the 
individual simulations computed with four bins would overall 
be relatively small. We performed test calculations with twelve 
opacity bins (with bin selection based on strength and wave- 
length), neglecting scattering extinction in optically thin layers 
or including it as true absorption, and found that the main results 
are very similar to the ones obtained with four bins. 

We also tested whether the cool stratification may be an arti- 
fact caused by the choice of mean temperature and density strat- 
ification in the definition of the mean opacities in the streaming 
limit (Eq. |5] and [TSl l. As mentioned in Sect. 12.31 we used for 
this task the mean stratifications from relaxed simulations that 
we computed excluding the scattering contribution to extinction 
in the streaming regime. A question is whether the computed 
streaming-regime bin opacities could force the temperature strat- 
ifications in the standard case to remain artificially cool - as well 
as in the case where we assume a Planckian source function and 
no scattering contribution to extinction in optically thin layers. 
To check this, we recomputed the streaming-regime bin opaci- 
ties for the metal-poor red giant at [Fe/H]=-3 for all three scat- 
tering treatments, instead using the hot mean stratification from 
the simulation in which scattering is treated as absorption ev- 
erywhere. We found consistency with the previous results: when 
scattering is correctly treated, or the contribution of scattering 
is left out from the streaming-regime opacities, the hot strati- 
fication immediately cools down to the same level that it had 
in the original simulation. In the scattering-as-absorption case, 
the temperature in the upper atmospheric layers of the simula- 
tions shows some degree of sensitivity to the mean stratification 
adopted for the opacity binning calculations. When the mean 
stratification for the opacity binning calculations is iteratively 
adjusted to match the actual mean stratification resulting from 
the scattering-as-absorption simulations, the temperatures in the 
surface layers eventually attain lower values than predicted in 
first approximation (e.g., about 70 K lower for the [Fe/H]=-3 
red giant simulation and about 200 K lower for the [Fe/H]=-3 
turnoff star), although they still remain significantly higher than 
in the coherent scattering and no-scattering-in-streaming-regime 
cases. Thus, our main conclusions are unaltered. We again stress 
that the scattering-as-absorption approximation is generally un- 
physical, and that the no-scattering-in-streaming-regime approx- 
imation provides a better approximation to the correct coherent 
scattering solution. 

Finally, Fig. [TO] shows the temperature-pressure distribu- 
tions in two representative snapshots from the [Fe/H]=-3 and 
[Fe/H]=-2 red giant simulations with proper treatment of coher- 
ent scattering in the radiative transfer. Although the distributions 
are limited from below by the curves defined by the process of 
H2 formation under adiabatic conditions, the stratification itself 
in the upper photosphere is actually not adiabatic. In all columns 
in the upper part of the simulations, the entropy gradient is gen- 
erally directed outwards - as expected for convectively stable 
regions - with only a few localized sites where the gradients are 
shallow or flat. Hence, even in the coherent scattering case (and 
also where the contribution of scattering to extinction in opti- 
cally thin layers is excluded), radiative transfer is still responsi- 
ble for regulating the temperature stratification in the uppermost 
layers, which would otherwise be even cooler if they were solely 
controlled by adiabatic processes. Consequently, we are confi- 
dent that the goodness of our approximation of radiative transfer 
with coherent scattering is not a spurious coincidence because 



the photospheric stratifications are possibly too close to the adi- 
abatic limit. 

5. Concluding remarks 

We have carried out 3D radiative hydrodynamic simulations 
of convection at the surface of three metal-poor stars (two gi- 
ants and one turnoff star) using different treatments of contin- 
uum scattering for the solution of the radiative transfer equation. 
These consist of two approximated implementations in which 
we assume a Planckian source function, excluding or includ- 
ing the contribution of scattering to opacity in the optically thin 
layers of the photosphere, and one implementation in which 
we self-consistently solve the radiative transfer equation for a 
non-Planckian source function that includes a coherent isotropic 
continuum scattering term. We found that simulations computed 
with the last-mentioned method and with the approximation in 
which we exclude scattering from the total extinction in opti- 
cally thin layers produce temperature and density stratifications 
that agree well with each other. This is essentially because, to 
first approximation, coherent isotropic scattering processes do 
not contribute directly to radiative energy exchanges in the up- 
per photosphere. In particular, both solutions predict a tempera- 
ture stratification in the upper photospheric layers of metal-poor 
late-type stars that is lower on average than the one predicted by 
ID hydrostatic model atmospheres. This result is consistent with 
the findings bv lAsplund et al.l (Il999h and lCollet etaP (l2007h : we 
therefore rule out that the cool temperature stratification pre- 
dicted by those works can be ascribed to their approximated 
treatment of scattering. 

The advantage of using the approximation of a Planckian 
source function and no continuum scattering contribution to ex- 
tinction in optically thin layers is that it allows us to achieve a 
sufficiently accurate representation of the radiative heating rates 
in the upper photosphere and produces temperature density strat- 
ifications very similar to the one obtained with the proper, iter- 
ative, solution, but at a much lower computational expense. In 
our analysis, we also showed that assuming a Planckian source 
function and treating scattering as true absorption everywhere in 
the simulation - as in present CO^BOLD models - leads to sys- 
tematically hotter temperature stratifications. At metallicities of 
[Fe/H]=-3 and [Fe/H]=-2, this approximation can produce de- 
viations from the correct solution of the order of a few hundred 
kelvins in the upper photosphere, which can significantly affect 
the predictions from spectral-line-formation calculations. 
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